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A jammed packing of frictionless spheres at zero temperature is perfectly specified by the network 
of contact forces from which mechanical properties can be derived. However, we can alternatively 
consider a packing as a geometric structure, characterized by a Voronoi tessellation which encodes 
the local environment around each particle. We find that this local environment characterizes 
systems both above and below jamming and changes markedly at the transition. A variety of 
order parameters derived from this tessellation carry signatures of the jamming transition, complete 
with scaling exponents. Furthermore, we define a real space geometric correlation function which 
also displays a signature of jamming. Taken together, these results demonstrate the validity and 
usefulness of a purely geometric approach to jamming. 


I. INTRODUCTION 

Over the past two decades the jamming of athermal 
frictionless spheres has been seen as the limiting case 
of several different kinds of systems [THS]. Athermal soft 
sphere systems can be brought to the limit of zero inter¬ 
nal energy and isostaticity, achieving a critically jammed 
system which is typically characterized by mechanical 
properties [a 0 imq]. However, when such systems 
are below the jamming density there is no longer a me¬ 
chanical network of force-bearing contacts and so me¬ 
chanical order parameters are all identically zero. Con¬ 
versely, hard sphere thermal liquids are studied below the 
glass or jamming transition and are characterized by dy¬ 
namic quantities such as mobility and pressure main]. 
As density is increased they reach the limit of diverg¬ 
ing reduced pressure and become a critically jammed 
system. Above this density, hard sphere systems can 
not exist. While both athermal soft sphere systems and 
thermal hard sphere glass systems have been successful 
models for predicting and measuring scaling exponents 
of various parameters near the jamming phase transition 
[aia Ellin], neither of these model systems speak to the 
behavior of unjammed athermal systems. This leaves a 
gap in the understanding of the athermal jamming transi¬ 
tion. In this paper we introduce new geometric order pa¬ 
rameters which characterize the athermal jamming tran¬ 
sition both above and below jamming, placing both sides 
of the transition on equal footing and providing a mean¬ 
ingful way to interrogate soft sphere systems below the 
jamming transition. 

The structure of jammed systems has long been stud¬ 
ied in terms of geometry HMD however a systematic 
study of geometric changes as a function of distance to 
the transition has not yet been performed. The Voronoi 
tessellation [TS] , which is well defined at all packing frac¬ 
tions, provides a natural lens through which to study 
both unjammed and overjammed systems. In previous 
work we have demonstrated that the number of facets 
(corresponding to the number of nearest neighbors) pro¬ 


vides a good order parameter for the jamming transition 
m- This order parameter raised a new problem, how¬ 
ever, because it exhibited an upper critical dimension 
(above which, all order parameters share the same scal¬ 
ing laws) of d = 3. This stood in contrast to the well 
known fact that mechanical order parameters exhibit an 
upper critical dimension of d = 2 [2110]. This, coupled 
with the recent success of replica theory in predicting 
high finite dimensional scaling |0] has motivated us to 
explore a range of geometric order parameters in dimen¬ 
sions ranging from d = 2 to d = 5. 

In this paper we show that most geometric properties 
of the Voronoi tessellation are controlled by the jamming 
point 4>j, suggesting that jamming can be described in 
purely geometric terms. Further, we present a new ge¬ 
ometrically defined correlation function which changes 
qualitatively at the jamming transition. Surprisingly, 
none of these measures show any indication of the previ¬ 
ously discovered pre-jamming transition, associated with 
the maximum inscribed sphere of the Voronoi cell, which 
we have found to happen at a density (p* < pj [El- 

II. GENERATING A PACKING 

We simulate packings of frictionless athermal particles 
with a harmonic contact potential in periodic boundary 
conditions as described in references [T9ll21]. In d = 3—5, 
we use monodisperse spheres, and in two dimensional sys¬ 
tems, we use a 50:50 mixture of bidisperse disks with a 
ratio of radii that is 1:1.4, known to show mechanical 
jamming. We present data obtained with three pack¬ 
ing protocols: Infinite Quench (IQ) (2, Geometric Mean 
(GM)[T2[1I], and Energy Sweep (ES)E]. 

Our three protocols differ only in how jamming is ap¬ 
proached. All begin with a set of particles in random po¬ 
sitions at a specified density. The energy of this system is 
then minimized to hnd the so-called inherent structure, 
found at the local energy minimum. Each of these proto¬ 
cols works as an iterative process by finding the inherent 
structure at a given density and then using this packing 
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FIG. 1. Plots of scaled order parameters vs. the scaled pack¬ 
ing fraction. Closed circles represent IQ data, x’s represent 
GM data (from below), and triangles represent ES data (from 
above). The parameters shown are (a) mean surface area, S, 
(b) standard deviation of volume divided by the mean of the 
volume, V, (c) mean surface to volume ratio (d) 

mean aspect ratio. A, and (e) mean aspect ratio angle cos9. 
We plot data for d = 2 (smaller particles light gray, larger 
particles dark gray, combined black), d = 3 (green), d = 4 
(red), and d = 5 (blue). 


as the seed to find a minimized packing at a new density. 

The IQ protocol begins with a random packing at 
zero density. Every particle is inflated to achieve a new 
packing at a specified higher density and this packings 
energy is then minimized. The results of this minimiza¬ 
tion are then used to create a denser packing, and so on. 
This proceeds in linearly spaced steps of packing fraction 
until the desired range of packing densities is covered. 
The range is chosen to cover densities from ^ = 0 to 
(j) = The limits of this range are somewhat arbi¬ 

trary but are chosen to be symmetric about (j)j. While 
the most relevant region is near the transition point, we 
include data at both the high and low extremes for com¬ 
pleteness. Data for d = 3 — 5 uses 65536 (2^®) particles, 
while d = 2 uses 16384 (2^"^) particles. 

The GM protocol is designed to zero in on the transi¬ 
tion point, either approaching from above or below, with¬ 
out ever overshooting. In this manuscript, we only report 
on GM systems approaching from below because the ES 
protocol (described below) converges much faster when 
approaching from above. The GM protocol requires an 
initial bounding of the jamming point by choosing two 
densities, one above and one below. A packing is ini¬ 
tially created at the lower bound and its energy mini¬ 
mized. A new packing is created between the upper and 
lower bounds using the original packing as its seed. If this 
packing is below jamming (taken to mean an energy per 
particle of < 10“^°), it becomes the new lower bound and 
serves as the next seed. If, however, this packing is found 
to be above jamming it is discarded and its density is used 
as the upper bound in picking a new intermediate den¬ 
sity. This proceeds until we approach the jamming point 
to within our energy per particle tolerance of 10“^°. In 
this way we are able to create a packing right at the edge 
of jamming that is the result of only inflationary steps, 
without ever crossing into the jammed regime. Because 
the convergence is slow, we are only able to report on 
8192 (2^^) particles. 

The ES protocol is limited in that it can only serve 
to approach jamming from above, but as previously men¬ 
tioned, it converges faster than GM. The ES protocol 
exploits the scaling of system energy with excess pack¬ 
ing fraction E (X {(j) — (j>j)^ to gently approach jamming 
from above, creating Usteps logarithmically spaced pack¬ 
ings per decade. Given an initial system density 4>i, sys¬ 
tem energy Ei, and a guess for the jamming density 4>i 
we calculate the packing fraction for the next system as 

= ( 1 ) 


Once this new system’s energy is minimized we compute 
a better estimate for the true jamming density as 


4>i+l 


4^i+l ^i/^i—1 

1 ~ \/Ei/Ei-i 


( 2 ) 


This process continues until we achieve an energy per 
particle of 10“^°. 
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FIG. 2. Log-log plots of each scaled order parameter vs. the scaled packing fraction approaching jamming from below (left) 
and above (right). Closed circles represent IQ data, x’s represent GM data (from below), and triangles represent ES data (from 
above). The parameters shown are (a,f) mean surface area, S, (b,g) standard deviation of volume divided by the mean of the 
volume, V, (c,h) mean surface to volume ratio (d,i) meau aspect ratio, A, and (e,j) mean aspect ratio angle cos9. 

We plot data for d = 2 (smaller particles light gray, larger particles dark gray, combined black), d = 3 (green), d — 4 (red), and 
d = 5 (blue). 
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We choose the starting point of the approach to be 
approximately 2(f)j. It has been previously shown that 
the jamming density when approaching from above is 
dependent on the initial packing density for systems that 
start close to phij [5T]. We choose to start at such a high 
value of (j) to ensure that our results are independent of 
the starting density. 

All ES data sets use 16384 (2^^) particles. Data for 
d = 2, d = 3, and d = 4 are averaged over 10, 63, and 79 
systems respectively while data for d = 5 is taken from a 
single run. 


III. GEOMETRY OF THE VORONOI 
TESSELLATION 



FIG. 3. Illustration of the aspect ratio axes in two Voronoi 
cells. For each cell, the short axis is shown in green (short 
dashes) and the long axis is shown in orange (long dashes). 
The angle 6 between the two axes is defined to be the acute 
angle between the short and long axis. The angle a between 
two long axes of different cells is also shown. 


Given a packing created via any of our protocols and in 
any dimension we calculate the Voronoi tessellation using 
the algorithms described in [TH] and extract the associ¬ 
ated vertices using the Delaunay triangulation [22]. For 
the monodisperse packings we create in d = 3 — 5, this 
Voronoi tessellation is the standard Voronoi tessellation 
wherein the size of a cell is independent of the size of the 
particle. However, due to the bidispersity used in d = 2 
we use the radical Voronoi (or Laguerre) tessellation [TB| , 
which makes the boundaries between cells the bisecting 
plane between the particle edges. This preserves the con¬ 
vexity of each cell and thus provides a natural extension 
of the classical Voronoi cell. From each Voronoi cell, we 
extract all of our measurements. The number of facets of 
the Voronoi tessellation gives us 1) the number of near¬ 
est neighbors V; The vertices of the Voronoi cell allow 
us to calculate 2) the surface area S and 3) the volume 
V; The ratio between the largest and smallest possible 
distances between parallel planes kissing the cell defines 
4) the aspect ratio A; Finally, the dot product between 
the headless vectors defining the aspect ratio provides 5) 
the cosine of the cell’s internal angle d. 


Parameter y 

N 

S 

V 


A 

cos d 

Power 7 

0.7 

1.0 

0.75 

1.0 

0.75 

0.33 

Xj, d = 2, large 

6 

3.005 

0.0315 

3.738 

1.221 

0.567 

Xj, d — 2, small 

6 

2.275 

0.0478 

3.827 

1.306 

0.559 

Xj, d = 2, all 

6 

2.640 

0.0573 

3.782 

1.264 

0.563 

XJ,d=‘i 

14.29 

5.385 

0.0386 

5.386 

1.322 

0.429 

XJ, d. = 4: 

32.74 

6.874 

0.0366 

6.875 

1.379 

0.385 

XJ,d = 3 

74.62 

8.261 

0.0340 

8.262 

1.412 

0.350 


TABLE I. Scaling laws and critical values for all parameters 
Y, such that 2t2C7 )^. All critical values are unitless 

except for Sj which is reported as the unitless 

For d = 2, we report separately on yj values for the larger 

particles, the smaller particles, and the system as a whole. 

A. Volume and Surface Area 

Calculating volumes and surface areas is notationally 
complicated but conceptually simple to achieve by break¬ 
ing the cell into simplices. To find volumes and surface 
areas we exploit the fact that the d-dimensional volume of 
a d-simplex can be calculated from the generalized triple 
product of its vertices. The Delaunay triangulation of the 
surface of a Voronoi cell breaks down the surface of each 
facet k into a number of (d — 1) dimensional simplices 
labeled by the index m. There are d-vertices associa¬ 
tion with each simplex, which we denote as Vm,i where 
i ranges from 1 to d, and we denote the outward facing 
normal vector to a facet k as fik- From this, the surface 
area of each facet is calculated as the sum of the surface 
of all of its constituent simplices as 

^ ' [('^771,1 '^7n,d) A * * * A {Vm,d—1 | 

'^^2^ (d-1)! ’ 

(3) 

where A denotes the d-dimensional wedge product. 

The total surface area of a given Voronoi cell is then 
the sum of all facets 

S = J2Sk- (4) 

k 

By choosing an interior point of the cell r, we can sub¬ 
divide the volume of the cell into a number of d-simplices 
whose volumes sum to the volume of the cell as 

TA K^m.d - fO • [(Dn.l - Vm,d) A • • • A {Vm,d-1 “ Dn,d)] | 

m 

(5) 

For a given packing, the mean cell volume is just the 
simulation volume divided by the number of particles. 

The distribution of cell volumes, however, does change. 






















5 


and so we report on V, the ratio of the standard devi¬ 
ation of the volume distribution to the mean. We also 
report on the mean of the unitless surface to volume ratio 


B. Aspect Ratio and Internal Angle 6 

The ratio of surface area to volume S/V^ defines a 
simple notion of an aspect ratio, but one that is insen¬ 
sitive to the anisotropy of the cell. We define another 
aspect ratio, explicitly sensitive to anisotropy by looking 
at the ratio between the longest one dimensional span 
in a cell to the shortest one dimensional span of a cell 
(Figure [^. To calculate this aspect ratio we define the 
long axis £ as the maximum distance between any pair 
of vertices and the short axis s as the minimum of the 
set of maximum distances between each vertex and each 
facet. Given a set of vertices vl and introducing a point 
Pk on each facet fc, these definitions can be formalized as 


£= {£ I 11^1 = MaxijWvi - {;j||}, (6) 

and 

s = {s I ||s|| = Minfe(Maxi|hfc • {vi -pk)\)}. (7) 

The aspect ratio is then simply defined as 



( 8 ) 


We can further measure the skewness of a cell by defin¬ 
ing the angle between the long axis and the short axis as 


cosd = 



(9) 


where the absolute value is taken because these are head¬ 
less vectors. 


C. Correlation Function 


We can examine the interaction of each cell with its 
neighbors by defining a correlation function based on the 
angle between the axes of pairs of cells. When cells are 
packed together to fill space neighboring cells must share 
facets, potentially causing the axes to align. To charac¬ 
terize this we measure the cosine of the angle between 
two long axes £i and £j associated with particles i and j 
respectively (illustrated in Figure [^. Because the axes 
are headless vectors we must use the formalism of direc¬ 
tors, giving rise to the definition for the cosine as 


cos ttij 


\\i\m\' 


( 10 ) 





(I nterparticle D istance )* N partkie 



(Interparticle Dlstance)*Npartde 


FIG. 4. The normalized long axis correlation between Voronoi 
cells plotted as a function of distance between particles in 
d = 2 — 5. Correlations are shifted vertically to show the 
effect of changing (j) with a color scheme that fades from purple 
{4> = 0) to black {(j) = (f)j) to green (0 = 20j). Gray dashed 
lines show the line corresponding to completely uncorrelated 
axes, open circles denote minima and open stars represent 
secondary maxima. A black line has been drawn over each 
curve representing the Savitsky-Golay filter. Data obtained 
using the IQ protocol. 
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To compare systems in different dimensions, we must 
first calculate the expectation values of completely un¬ 
correlated directors. The expectation value of the cosine 
of the angle in dimension d is given by 


(cOSQ!)d = 


r7r/2 . d—2 i 

Jq cos a sin ada 


r(f) 


Sin 


' da 




( 11 ) 


The standard deviation of the angle between uncor¬ 
related directors in dimension d is defined as ad = 
yj (cos — (cos^ a)d- Therefore we also calculate the 
expectation of the square of the cosine of the angle of 
uncorrelated directors as 


(cos^ a)d = 


cos a sin 


' ada 


r 

Jo 


sm 


' da 


r(f) 

2r(d+i)- 

( 12 ) 


Thus we find the standard deviation of uncorrelated di¬ 
rectors in dimension d to be 


= 


r(f) 


r(i) 


\ 7rr(^^)" 2r(d+i)' 


(13) 


We define our correlation function as the normalized 
value of the cosine of the angle between the long axes 
of every pair of particles as a function of the distance 
between cells as 


Ce{r) 


1^ _ II ^ cos aij - {cos a)d 

''‘-oW-O - ’—-I -■ («) 


We note that this correlation function, relating the 
shape and asymmetry of Voronoi cells, is logically dis¬ 
tinct from the pair correlation function, or indeed from 
any correlation function based solely on particle posi¬ 
tions. 


IV. ANALYSIS 


A. Order Parameters 

Figure [2 presents the geometric order parameters de¬ 
scribed above calculated for systems created with all pro¬ 
tocols as a function of distance to (/>j for d = 2 — 5. Data 
across multiple dimensions is presented on the same scale 
by subtracting off the value at the jamming transition 
and then dividing by that same value. The packing frac¬ 
tion is similarly scaled as {(j)—(j)j)/4>j- Below jamming all 
of these parameters change rapidly with increasing pack¬ 
ing fraction. Jamming is marked by a sharp kink and 
above jamming they evolve with a much gentler slope. 
For all measures except the surface area to volume ratio 
the d = 2 data does not seem to collapse onto the same 
family of curves as the higher dimensions. Because the 


d = 2 packings are bidisperse we show separate curves 
for each particle size (shown in dark and light gray) and 
a single curve representing the combined data (shown in 
black). The difference is especially apparent in the sur¬ 
face area: the Voronoi surface area increases for larger 
particles and decreases for smaller particles as jamming 
is approached from below. This makes intuitive sense; 
at extremely low packing fractions the Voronoi cells for 
the two sets of particles should be almost identical and 
near jamming the larger particles will end up having a 
larger surface area and a larger volume than their smaller 
counterparts. In the combined data, the curve collapses 
to follow the trend observed in d = 3 — 5. 

In order to explore the behavior very close to jamming 
we use the GM protocol to approach from below and ES 
to approach from above. In this way we obtain packings 
that converge logarithmically on (j)j. Plotted on a log-log 
scale (Figure we find that each parameter scales with 
its own power-law on both sides of the transition with 
power law values and critical values listed in Table We 
have previously demonstrated that the mean number of 
neighbors (TV) scaling is consistent with a power of ~0.7 

m- 

Below jamming, there must be a limit to the scaling 
regime. The mean surface area (S), volume (V), and 
number of facets (TV) of Voronoi cells at ^ = 0 and 
their respective dimensional dependence can be semi- 
analytically determined O . The same should be 

true for aspect ratio (A) and internal angle (cos 6*) but 
to our knowledge those studies have not yet been done. 
This is responsible for the changes in curvature seen at 
low (p in Figure 

While most of the power laws work well over at least 
five decades, there are a few exceptions. The data from 
below is very sparse, and so we cannot claim that the 
power laws fit exactly and can only suggest that the plots 
look like power-laws within the plotting area. Precise 
claims about the scaling exponents of these power laws 
would require a method which approached jamming more 
predictably from below and which converged much faster 
so that averaging could be used, as it is done above jam¬ 
ming. 

It is also important to note that the d = 2 data devi¬ 
ates significantly in the standard deviation of the volume 
(Figure [^) and the internal angle (Figure [^). When 
coupled with the fact that the mean number of neigh¬ 
bors does not show a signature of jamming in d = 2 m, 
this strongly suggests that d = 2 is below the upper crit¬ 
ical dimension of the jamming transition when viewed 
from a geometric perspective. 


B. Correlation Function 

From the measurements of the aspect ratio we can 
see that at jamming the Voronoi cells are much more 
isotropic than they are far from jamming. At jamming, 
the aspect ratio is close to 1 and the direction of the long 
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FIG. 5. The position of the minimum (closed circles) and 
secondary maximum (open star) of the correlation shown in 
Figure as a function of distance from the jamming transi¬ 
tion. Colors shown represent dimensions 2 (black), 3 (green), 
4 (red), and 5 (blue). 

and short axes are uncorrelated as measured by cos 9. 
In contrast, at very low density the cells are elongated 
and have a large aspect ratio and axes that are nearly 
perpendicular. Figure]^ shows the measured correlation 
function between the long axes as a function of interpar¬ 
ticle distance for packing fractions ranging from ^ = 0 
to (j) = 2(^j in dimensions d = 2 — 5. Far below jam¬ 
ming, neighboring particles are highly correlated. This 
correlation decreases with increasing distance, showing 
an anti-correlated dip at intermediate distances and then 
finally decaying to completely decorrelated at large dis¬ 
tances. At jamming, the correlation function changes 
qualitatively, marked by the appearance of a positive cor¬ 
relation peak at intermediate distances in addition to the 
short distance dip. Both the dip and the peak become 
more prominent and sharpen at higher packing fractions. 
These extrema are found using a cubic Savitzky-Golay fil¬ 
ter with a span of 51 data points |26j and the positions 
of the dip and peak are indicated by circles and stars re¬ 
spectively in Figure]^ and plotted as a function of pack¬ 
ing fraction in Figure The position of the maximum 
shows a clear signature of the transition in d = 2 — 5. 
However, the position of the minimum for d = 3 — 5 


does not show a clear signature of this transition. We 
find that the correlation functions plotted in Figure 
only depend on interparticle distance, with no angular 
dependence when oriented to the long axes of each given 
particle. This correlation function is unusual in that the 
jamming transition is marked by the disappearance of 
the nearest neighbor correlation, seen in the value of the 
correlation function at the shortest possible interparticle 
distance. 


V. CONCLUSION 

We have observed a clear signature of the jamming 
transition in each of the studied measures of the Voronoi 
cell as well as in our newly defined axis-correlation func¬ 
tion. These results bolsters the claim that while jamming 
is a mechanical transition, it can be viewed separately as 
a purely geometric phenomenon. These results justify 
the use of the Voronoi cell as a tool to understand the 
jamming transition. Furthermore, we provide evidence 
that while the mechanical transition has an upper criti¬ 
cal dimension of d = 2, the geometric transition has an 
upper critical dimension of d = 3 when considering some 
geometric order parameters. 

Ultimately, each of the measures are sensitive to the 
fluctuations in the size and shape of individual Voronoi 
cells. Each measure reflects a different change in the 
cell. The fact that we see power-law scaling in all of 
these measurements suggests that nearly every aspect of 
the cell changes and is controlled by the transition from 
unjammed to jammed. Our results demonstrate that the 
mechanical jamming transition coincides perfectly with 
a transition in the geometry of the packing at (fij. 
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